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Abstract 

We discuss a simple extension of the Ho and Lee model with generic time- 
dependent drift in which: 1) we compute bond prices analytically; 2) the yield 
curve is sensible and the asymptotic yield is positive; and 3) our analytical 
solution provides a clean and simple way of separating volatility from the 
drift in the short-rate process. Our extension amounts to introducing one 
or two reflecting barriers for the underlying Brownian motion (as opposed to 
the short-rate), which allows to have more realistic time-dependent drift (as 
opposed to constant drift). In our model the spectrum - or, roughly, the set 
of short-rate values contributing to bond and other claim prices - is discrete 
and positive. We discuss how to calibrate our model using empirical yield 
data by fitting three parameters and then read off the time-dependent drift. 
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1 Introduction 


An allure of short-rate models is their apparent elegance and simplicity. In theory, 
once the short-rate process and the corresponding risk neutral measure are speci- 
hed, bonds and other claims, such as bond options, can be priced via simple-looking 
conditional expectations of “exponentially” discounted claims. However, in practice, 
things are a bit trickier. In their simplest incarnations, e.g., Merton’s (1973) model 
and the Ho and Lee (1986) model, short-rate models typicalljl^] produce unrealistic 
yield curves for zero-coupon bonds, with long-maturity yields turning negative. 

An immediately evident - but not necessarily relevant - culprit would appear to 
be that in such models can become negative. However, intuitively, since is not 
observable in real life, it taking negative values need not be a big deal. Instead, what 
is important is that the model yield curve be sensible. Thus, in the Vasicek (1977) 
model Tt can become negative, but for a range of parameters the model yield curve a 
priori appears sensible as the spectrum in this model is discrete and bounded from 
below. In models such as Merton’s model and the Ho and Lee model the spectrum 
is continuous and unbounded. Intuitively, the “spectrum” here can be thought of as 
the set of values of r* that actually contribute to the bond (and other claim) prices0 
While the spectrum in the Vasicek model is bounded from below, it is nonnegative 
only if, roughly speaking, volatility is low. In some cases, including in the current 
low interest rate environment with not-so-low volatility, this can be problematic. In 
fact, one can ensure that rt is always positive and can never reach zero as, e.g., in the 
Black and Karasinski (1991) model, where ln(r 4 ) is the Vasicek process. However, 
here too, low interest rates by construction imply low (almost vanishing) volatility. 

An alternative approach to ensure ri > 0 is to introduce a reflecting barrier at 
rt = 0, or to treat rt as an option on an underlying “shadow” rate Xt (which can be 
negative) by taking its positive partj^as in (Black, 1995) and (Rogers, 1995, 1996). 
Short-rate models with a reflecting barrier at r^ = 0 were studied in (Goldstein 
and Keirstead, 1997) for Merton’s model, the Vasicek model, the extended (a.k.a. 
shifted) CIR (Cox, Ingersoll and Ross, 1985) model, and the Longstaff (1989) model. 
Models with r* as an option were studied in (Gorovoi and Linetsky, 2004) when the 
“shadow” rate Xt follows the Vasicek model and the shifted CIR model. In this 
approach volatility generally need not be small even in the low interest rate regime. 
In all the models studied in (Goldstein and Keirstead, 1997) and (Gorovoi and 

^ One can circumvent this via apparently unrealistic time dependence for the drift in rj. 

Mathematically, the “spectrum” is the set of eigenvalues for the static Schrodinger equation 
to which the pricing PDE reduces, e.g., when the underlying parameters are time-independent. 

® For a reflecting barrier at = 0 we have a Neumann boundary condition w.r.t. rt {rt G R^) 
and the system can be described via rt = \Xt\ with claims symmetric under Xt —>■ —Xt. The 
spectrum is discrete and positive even if W € R follows Merton’s model or the Ho and Lee model 
(see below). For r* as an option we have rt = (2ft)''' (with the “shadow” rate Xt G R), and at 
rt = 0 we sew together solutions with rt = 0 (Xt G R” ) and rt = Xt G R''" by requiring continuity 
of the claim price and its first derivative w.r.t. Xt. Here the rt = Xt G R model is assumed to 
have a discrete spectrum, and then the spectrum in the rt = (2ft)''' model is discrete and positive. 
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Linetsky, 2004) the underlying parameters (volatility, drift, mean-reversion rate, 
etc.) are assumed to be constant. More generally, in time-homogeneous cases (where 
the parameters depend on (or Xt) but have no explicit time dependence), we 
have the standard separation of variables and the pricing PDE reduces to an ODE 
(static Schrodinger equation). Also, at = 0 we have a time-independent boundary 
condition for a reflective barrier and a time-independent sewing condition for as 
an option. These nice simplifying features are lost once we consider time-dependent 
parameters, which in many cases are required to describe real-life yield curves. 

However, not all is lost. In this paper we set forth a simple way of circumventing 
this difficulty. The following observations will pave the way for us. First, to get a 
sensible yield curve (and, more generally, claim pricing), it is not required that rt be 
nonnegativel^ What is required is that the spectrum be discrete and nonnegative0 
To achieve this, it suffices to introduce a reflecting barrier for the underlying Brow¬ 
nian motion Wt (as opposed to r*). Then we can have time-dependent drift and 
still achieve separation of variables. Furthermore, the resulting boundary condition 
is also time-independent. Second, with a reflecting barrier for Wt, the spectrum is 
discrete and positive (roughly, if rt > r_ < oj^ already in the Ho and Lee mode@ 

drt = a dWt + vit) dt ( 1 ) 

where we take constant volatilit y a and general time-dependent drift z/(f) and solve 
the pricing problem analytically!^ Because of the simplicity of the diffusion part in 
the Ho and Lee model, our analytical solution provides a clean and simple way of 
separating volatility in rt from the drift. We discuss how to ht this model into an 
empirical yield curve, which involves calibrating three constant parameters (volatil¬ 
ity cr, the initial value of rt, and the location of the reflecting barrier) and then read 
off the time-dependent drift, or, equivalently, its contribution to the yield curve. We 
illustrate our method by htting the model into recent U.S. Treasury yield data. 

The remainder of this paper is organized as follows. In Section 2 we briefly review 
some generalities of short-rate models and then discuss how to introduce reflecting 
barriers for short-rate processes of the form rt = U{Wt) + x(t) and solve the claim 
pricing problem. In Section 3 we apply the results of Section 2 to the Ho and Lee 
model with one and two reflecting barriers, give explicit formulas for zero-coupon 
bond prices, and for illustrative and comparison purposes discuss htting the model 
into the Japanese Government Bond data used to calibrate a model in (Gorovoi 
and Linetsky, 2004). We also discuss a ht into recent U.S. Treasury yield data. We 
briehy conclude in Section 4. Some technical details are relegated to Appendices. 

® Thus, the models studied in (Goldstein and Keirstead, 1997) and (Gorovoi and Linetsky, 
2004) can have: a reflecting barrier at = r_ < 0; and rj = r_ -|- (W)'*’ with r_ < 0, respectively. 

^ Realistically, the lowest eigenvalue has to be reasonably positive (see below). 

® I.e., for this purpose alone we do not need more complex (Vasicek, GIR, etc.) dynamics. 

® More generally, we discuss factorized processes of the form rt = U{Wt) + x{t)- Furthermore, 
if desired, we can have ct > 0 by appropriately choosing parameters. However, this is not required. 

Both for one reflecting barrier (say, at Wt = 0) and two reflecting barriers (at Wt = 0 and 
Wt = L). The latter case can be useful for preventing rt from wandering away into large values. 
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2 Short-rate Models 


A short-rate model posits a risk-neutral measure Q and a short-rate process r^. The 
cash bond process is given by 


Bt 



( 2 ) 


while the price at time t of a claim X at maturity T is given by v{rt,t,T), where 
the pricing function v{z,t,T) is given by 

v{z,t,T) =/exp (- f Ts ds\ x\ (3) 

\ \ Jt J I (^^rt=z 

E.g., for X = 1, we have v{rt,t,T) = P{t,T), and v{z,T,T) = 1, where P{t,T) is 
the price of a zero-coupon T-bond. Also, {■) denotes expectation. 

In short-rate models one usually works with a parameterized family of processes, 
and chooses the parameters to best £t the market. Thus, typically one assumes that 
Tt satishes the following SDE: 


dvt = a{rt, t)dWt + virt, t)dt (4) 

where (j{z, t) and h>{z, t) are deterministic functions, and Ifd is a Q-Brownian motion. 
The pricing function v{z,t,T) satishes a pricing PDB^T which follows from the 
requirement that the discounted process Z{t,T) = B^^v{rt,t,T) be a martingale 
under the risk-neutral measure Q; 

dtv{rt,t,T) + - rtv{rt,t,T) = 0 (5) 

with the terminal condition v{z,T,T) = {X)ci^rT=z = 


2.1 Mean-Reversion and Positivity 

In some cases one may wish to require that the short-rate process not wander 
away to large values. One way to achieve this is to use a mean-reverting process 

dXt = a{t)dWt + [6{t) — a{t)rt] dt ( 6 ) 

where cr(t), 6{t) and a{t) depend only on time. For constant a, 6 and a we have the 
mean-reverting Ornstein-Uhlenbeck process. In the Vasicek/Hull-White model we 
have Tt = Xt- One “shortcoming” of this model is that rt can occasionally become 
negative^ One way to deal with this is to consider short-rate models of the form 
A = ■'"o /(-^t)//(0), where f{y) is a positive function, e.g., f{y) = exp{y), which 
is the Black-Karasinski model. The path integral treatment of such models was 
discussed in (Kakushadze, 2014) @ 

Known as the Feynman-Kac equation, a.k.a. the Arrow-Debreu security PDE. 

More relevantly (see below), e.g., for constant a, 6 and a, the asymptotic {i.e., large T) 
zero-coupon bond yield in this model is positive only for P < 2a9. 

The results we obtain below can also be derived using path integral. 
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2.2 Short-rate Models with Reflecting Barriers 

An alternative approach is to consider short-rate processes of the form rt = UiWt-, t), 
where the function U{x, t) is bounded. More precisely, here we can consider three 
cases: i) x is unrestricted, a; G R; ii) x takes values on a half-infinite line, x G [a;_, cxd) 
or X G (— cxD,x+]; and hi) x takes values on a hnite interval, x G [x_,x+]. In the 
case hi), U{x,t) can be a simple function, e.g., U{x,t) = a{t)x + b{t), even though 
it is not bounded when extended to the entire real line x G R. In the case ii), r+ is 
nonnegativS so long as U{x, t) > 0 for, say, x G [x_, oo). 

In the following we will assume that, for the allowed values of x, for a given 
value of t, there is a one-to-one mapping between z and x via z = U{x,t). Then the 
pricing PDE (jS]) simplifies as follows. Let Y{U{x, T)) = Y (x) ancii^x(17(x, t),t, T) = 
'ip{x,t). Also, we have a{rt,t) = dxU{x,t) and v{rt,t) = d‘lU{x,t)/2 + dtU{x,t), 
where U{x,t) = rt, and the pricing PDE ([5]) reads: 

dti){x,t) = ~ dl'tp{x,t) + U{x,t) 'ip{x,t) (7) 

ij{x,T)=Y{x) (8) 

where ([8]) is the terminal condition at t = T. For definiteness, here we focus on 
the case with two boundaries. Cases with a single (lower or upper) boundary follow 
upon removing the unwanted boundary to the corresponding infinity and requiring 
that 'ip{x,t) be hnite (actually, vanish) at such inhnity. 

Since we have boundaries, we need to specify boundary conditions. For the 
reasons which will become evident momentarily, we will take the boundaries x± to 
be rehecting barriers, he., when the Brownian motion touches the lower (upper) 
boundary x_ (x+) from above (below), it is rehected back into the values x > x_ 


(x < x+). This implies that we have Neumann boundary conditions 

d^'4){x-,t) = t) (9) 

9^'0(x+,t) = O (10) 

This then implies that Y (x) also satishes Neumann boundary conditions: 

d^Y{x_) = 0 (11) 

d,Y{x+) = 0 ( 12 ) 

Had we imposed the Dirichlet boundary conditions 'ip{x-,t) = 0 and 'ip{x+,t) = 0, 


for which the process Wt is not allowed to touch the boundaries, we would invariably 
have F(x_) = 0 and E(x+) = 0. This would not be suitable for our purposes here 

Albeit, as discussed above and below, nonnegativity of rt is actually no longer required in 
cases ii) and iii) and is replaced by a weaker condition. 

When U{x,t) has no explicit time dependence, we have time translational invariance and 
'i/’(x,t) depends on T only via the combination T — t. However, this is no longer the case when 
dtU{x,t) ^ 0. Nonetheless, below we will use the abbreviated notation tj}{x,t). 
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as Dirichlet boundary conditions are incompatible with, e.g., the claim X = 1 for a 
zero-coupon T-bondO 

Note that ([7]) is the Schrodinger equation in imaginary (Euclidean) time for a 
particle (with mass m = 1 and Planck’s constant h = 1) in the potential U{x,t). 
For general time-dependent potentials U{x,t) Eq. ([7]) is difficult to solve. Howeyer, 
for our purposes here it will suffice to consider factorized potentials of the form 


U{x,t) = V{x)+x{t) (13) 

For such potentials we have = V'{x) and p{rt,t) = V''{x)/2 + x(t), where a 

prime denotes a derivative w.r.t. x, while a dot stands for a derivative w.r.t. t, and 
n = V{x) + x{t)- 

For the factorized potential flT^ . we have separation of variables, and the solution 
to (17|) can be written in terms of a series {Ei < F 2 < F 3 < ... )0 


00 



n=l 

( 14 ) 


vit,T) = j ds [x{s) - x{t)] 

( 15 ) 


Xn{t) =x{t)+ En 

( 16 ) 

where iJnix), n G N 

is the complete orthonormal set of solutions to the static 

Schrodinger equation: 

dl'tpn{x) +V{x) 'ijjnix) = En 'ijjnix) 

( 17 ) 


rx+ 

/ dx ljjn{x) iJn'ix) = Snn> 

J X- 

( 18 ) 

subject to the Neumann boundary conditions 




( 19 ) 


= 0 

( 20 ) 


If we have a single lower (upper) reflecting boundary at X- (a;+), then we have only 
one Neumann boundary condition at this boundary together with the requirement 
that ilJni.x) vanish as a; — )■ -|-cxo (—cxo). 

For a zero-coupon T-bond we have X = 1 and u(rt, t, T) = P{t, T), where P{t, T) 
is the bond price. The yield is given by 

T-t 

We could consider inhomogeneous boundary conditions such that Y{x) would not need to 
satisfy Neumann boundary conditions. However, this will not be needed for our purposes here. 

Here we are assuming that the eigenvalue spectrum En is discrete and bounded from below. 
This is the case if the potential V(x) is confining, even without the boundaries. With one or two 
boundaries, it suffices that V(x) is bounded from below for the allowed range of x. 
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The asymptotic yield must be nonnegativJ^ 


0 < = lim R{t, T) = x*{t) + xi{t) = x*{t) + x(t) + Ei 

T—t^+oo 


where 


X* 


lim 

T—t—y-\-oo 


’l{t,T) 

T-t 


Assuming x*if) is hnite, we have a requirement that Ei > —x*(t) — xif)- 
The coefficients c„ (for a general claim are given by 


nx+ 

Cn= dx ‘Ipnix) Y{x) 

J X- 


( 22 ) 

(23) 


(24) 


which is a consequence of ([H]) and flTHD . To recap, the pricing function v{z,t,T) is 
given by 

OO 

v{z,t,T) = (25) 

n=l 

where x is related to z via V (x) + x(t) = z. 


3 Ho and Lee Model with Reflecting Barriers 


Let us consider the simplest example, the Ho and Lee model 


dvt = a dWt + z/(f) dt 


(26) 


where a is constant, but the drift i>{t) a priori is an arbitrary function of t (subject 
to some restrictions we discuss below). We have 

rt = aWt + x(t) (27) 


where (we are assuming ITo = 0) 


X{t) = ro + 


ds h'(s) 


So, V{x) = a X. 

The static Schrodinger equation (1T7)1 reads: 


(28) 


- 2 + a X i)n{x) = En iln{x) 


(29) 


The solution depends on the boundary conditions and is expressed via a linear 
combination of Ai {ax — e„) and Bi {ax — e„), where a = (2cr)^/^, = E^/fd, (3 = 
a/a = ((j^/2)^/^, and the Airy functions Ai{y) and Bi{y) are the two independent 
solutions of the Airy equation f"{y) = y f{y)- 

Realistically, the asymptotic yield should not be less than some positive number. 
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3.1 Ho and Lee Model on a Semi-infinite Line 

Thus, let us consider the Ho and Lee model fl26l) on a semi-inhnite line, he., we 
restrict the values of Wt to R'*', where without loss of generality we have set X- = 0. 
We have fl2^ with only one Neumann boundary condition 

d,A(0) = 0 (30) 

and the requirement that iJnix) be finite as a; —)■ +cxo. The solution is given by: 

iJuix) = an Ai{a X - e^) (31) 


where 

en = (32) 

Here (0 > > ^2 > • • •) are the zeros of the first derivative of Ai{y): 

At'iU) = 0 (33) 

Note that the integration in fl2T)) is from 0 to + 00 . The normalization coefficients 




dx [Ai {ax - en)f = -a ^^n [Ai (^n)]' 


(34) 


where we have used flTOj) (see Appendix IXj) . 


3.1.1 “Modulus” Ho and Lee Model 

In the Ho and Lee model on a semi-inhnite line the short-rate process rt generally 
is unbounded from above. If this is not problematic!^ then we do not even need 
to restrict the process to a semi-inhnite line. Instead, we can simply consider a 


“modulus” modelS 

n = a \Wt\ + x{i) (35) 

where Wt is unrestricted. In this case we have 

'ipn{x) = ttn Ai {a X — Cn) , o: > 0 (36) 

^n(x) = (-l)"+Vn(-a:), X < 0 (37) 

poo 

= 2 / dx [Ai {ax — e„)]^ (38) 

Jo 

a~'^ = 2 a~^en[Ai{-en)f , n = l,3,... (39) 

a-^ = 2 a-^ [Ai'{-en)f , n = 2,4,... (40) 


E.g., for given volatility cr, we are interested in time horizons such that rt simply does not 
have enough time to wander away too far. 

Note the difference with, e.g., an alternative model rt = IcrWt -|-x(t)|. The latter is harder 
to tackle for general x(t). See Appendix iBl for details. 
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and Bn are now given by 


( 41 ) 

( 42 ) 


Cn 'b(n+l)/2) 1,3,... 

B-n ini’ll ^ 2, 4, . . . 

Here Cn (0 > Ci > C 2 > • • •) are the zeros of 

Ai{U) = 0 (43) 

Note that the integration in fl241) is now from —00 to +cxd. 

The solutions 'ipn{x) with odd n are symmetric w.r.t. the Z 2 reflections x —)■ 
—X and satisfy the Neumann boundary condition fl30l) . while the solutions with 
even n are antisymmetric and satisfy the Dirichlet boundary condition V’n(O) = 0. 
The short-rate process is invariant under Wt —Wt, so the claims should also be 
symmetric: Y{—x) =Y{x). Then we have c„ = 0 for even n, and for odd n we have 

poo 

Bn = 2 dxipn{x)Y{x) (44) 

Jo 

It then follows (taking into account the differing by \/2 normalizations for ijjnix) 
in the two cases) that, for the same (symmetric) claim X = Y{Wt), the pricing 
function v{z, t, T) is the same in the Ho and Lee Model on a semi-inhnite line and in 
the “modulus” model fl35D - as they should be based on symmetry considerations. 


3.2 Ho and Lee Model on an Interval 


Assuming the reflecting barriers at a; = 0 and x = L, the Brownian motion Wt 
wanders between 0 and L. We have fl2^ with two Neumann boundary conditions 

dxi>ni0) = d^i>niL) = 0 (45) 

The solution to is given by 


^jJn{x) = ttn (t)n{x) 

(f)n{x) = Ai {a X — 

_ Ai' {a X 
Q{x,b} = 


Bi' (a X 


Cn) - <5(0, e„) Bi (cK X 
-e) 

-e) 


where e„ are the roots of the following equation for e: 


^n) 


( 46 ) 

( 47 ) 

( 48 ) 


Q{L,e) = Q{0,e) 


( 49 ) 


which has an inhnite tower of discrete solutions e„ > 0, n G N. The normalization 
coefficients 

= [ dx (pl{x) = - [Bn 0^(0) + {a L- Bn) 0^(L)] (50) 

io « 

where we have used fl4^ and (170]) (see Appendix]^. Note that if we take L —)■ cx), 
fl50|) reduces to (1341) . as it should. 



3.2.1 “Periodic” Ho and Lee Model 

Just as in the case of the Ho and Lee model on a semi-inhnite line, the Ho and Lee 
model on an interval is also equivalent to a model where Wt is unrestricted. In this 
unbounded model the short-rate process has the following “periodic” form 


rt = a h{Wt) -f xif) (51) 

where h{x) is a piecewise linear periodic function 

h{x) = X, 0 < X < L (52) 

h{x) = 2L — X, L < X <2L (53) 

h{x + 2L) = h{x) (54) 


and the claim X = Y {Wt) is symmetric under both x —)■ —x and L + x ^ L — x 
reflections: Y{—x) = Y{x) and Y(L — x) = Y{L + x). This implies that F(x) is 
periodic: Y{x + 2L) = Y{x). Therefore, the integration in (|2TD is from 0 to 2L, or, 
equivalently, from —L to L. 


3.3 Bond Pricing on a Semi-infinite Line 


For illustrative purposes, let us discuss zero-coupon T-Bond pricing, for which the 
claim is simply X = 1. For the Ho and Lee model on a semi-inhnite line we have 

noo nco ^ 

Cn= j dx iJn{x) = - an dy Ai (y) = - an Ai{^n) Gi'{^n) (55) 

Jo ^ Ji-a “ 

where we have used flTTl) (see Appendix E]) , and Giijj) is a Scorer function. 

The T-bond pricing function therefore is given by the following simple formula: 


v(^z,t,T) 




OO 


n=l 


TT Gz'i^n) 

I'Cnl Ai(^^n) 


Ai 



g-Xn(t)('r-b 


(56) 


where 

Xnit) = Xit) + En = Xit) + fJ \^n\ (57) 

and ri{t,T) is dehned in ffT5ll . Recall that /3 = (cr^/ 2 )^/^. 

Taking into account the asymptotics given in Appendix El we see that the series 
is well-behaved at large n. Thus, let 


v(^z, t,T) = Y^ Vn(z, t, T) (58) 

n=l 


For large n such that |^n| 3> ( 2 : — x(t)) //d, the leading asymptotics are: 

Vniz, t, T) ~ Un exp (-7 (t, T) 





2/3 


(59) 

(60) 
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Y 


(61) 




So, asymptotically, we have an alternating series, which converges according to the 
Leibniz criteriono In numerical computations one would truncate the series at a 
suitably chosen hnite n (see below). 

3.3.1 Parameter Count and Model Calibration 

Suppose we have data for zero-couporl^ T-bond prices. How many parameters would 
we need to £t to calibrate fl56|l ? We need to break this down. Let us start with the 
case of vanishing drift u{t) = 0. Then the answer is that we have 3 parameters to 
fit. Indeed, in fISB]) we have z, (3 (which is fixed by a) and tq (note that x(t) = tq 
when iy(t) = 0). However, looking at ([3]) for this model with iy(t) = 0 

dtv(z, t,T) + ^ cr^ t,T) - z v{z, t,T) =0 (62) 

v{z,T,T) = l (63) 

it might appear that we have only 2 parameters, and a, to £t. The third parameter, 
ro, is simply the location of the reflecting boundary^ which is a freS parameter. 
In fact, in the Ho and Lee model, before introducing a reflecting boundary, the 
spectrum of (162]) is continuous^ Once we introduce a reflecting boundary, the 
spectrum is discrete and, once again, its nonnegativity does not require that the 
reflecting boundary be at r_ = 0. It suffices to require that r_ = ro > —(3 

When the drift v{t) is nonzero, we haye a choice. Thus, we can try to £t more 
parameters. E.g., we can assume that z/(t) = z/q 7 ^ 0 is constant; then we have four 
parameters to £t. Or we can assume that u{t) is a general polynomial of degree 
k, so we haye 4 + k parameters to £t. Similarly, we can assume that v{t) is some 
function, e.g.^ v{t) = z/q cos(ci; t), in which case we haye hve parameters to £t. Etc. 

Alternatively, we can follow a different procedure, which we set forth here. In¬ 
stead of trying to £t z^(t), we can hrst £t the three parameters z, (3 and tq {e.g., via 
the least squares method - see below) assuming z/(t) = 0 , and then attribute the 
difference between the so-£tted model yield curve Rm{t,T) and the empirical yield 
curve Re{t,T) to nontrivial z/(t). We will refer to this difference as the “residual” 
yield: Rr{t,T) = Re{t,T) — Rm{t,T). In fact, looking at fl56l) it is evident that for 
nontrivial v{t) it is more convenient to £t 2 ; (which is nothing but r^), (3 and x(t) (as 
opposed to ro). So, below, after we £t y, (3 and tq in the z/(t) = 0 case, we will use 

Note that Un = v^2/3 (1 —\/2) C(l/2) ~ 0.494, where C(s) is the Riemann zeta function. 

In practice, unless the yield data is readily available, we may have data for coupon-bearing 
bond prices, which then are bootstrapped to obtain zero-coupon bond prices. 

When v{t) = 0, the reflecting boundary at a:_ = 0 for the Brownian motion process W* 
translates into a reflecting boundary at r_ = tq for the short-rate process r*. As discussed in more 
detail in Appendix [B] this is no longer the case for a general time-dependent drift v{t). 

Modulo the requirement (l22l) . that is (see below). Note that x*(t) = 0 when v{t) = 0. 

Moreover, (l62ll has a symmetry under the transformation v{z,t,T) —>■ e~‘’^ v{z + (^,t,T) = 
V(^{z,t,T), where ( is an arbitrary constant. I.e., if v{z,t,T) satisfies (l62ll and (|63ll . then so does 
VQ{z,t^T) for arbitrary C. 
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the so-obtained vq as x(i) iii the nontrivial h>{t) case@ Then we have the following 
simple formula for the “residual” yield: 

= ( 64 ) 

SO we can read off the empirical yield curve once we £t Rm{t,T), which will 

give us approximate^^! shapes of x(s) and u(s) for t < s < T. 

3.3.2 An Illustrative Example: Japanese Government Bonds 

For illustrative and comparison purposes, we have used the zero-coupon T-bond 
pricing formula in the Ho and Lee model on a semi-inhnite line with = 

0 (“Model-2”) to £t the Japanese Government Bond data used in (Gorovoi and 
Linetsky, 2004) to calibrate Black’s model of interest rates as options (Black, 1995)@ 
with the underlying “shadow rate” following the Vasicek model (“Model-1”). As in 
(Gorovoi and Linetsky, 2004), we £t the model by minimizing the root mean squared 
error (RMSE) between the empirical yield curve (Table (H column 4) and the Model-2 
yield curve (see above). 

Table [1] summarizes our results. The calibrated model parameters in Model-2 
ar^ 2 —0.00184, (3 ~ 0.0924 (which implies a ~ 0.0397), tq —0.05834. This 

implies that the initial short-rate value rt = z (here t = 02/03/2002) is essentially 0 
(within the 2-digit precision of the underlying bootstrapped yield data in column 4 of 
Tabled]). On the other hand, in Model-2 the lower bound on the short-rate process is 
r_ = ro ~ —5.834%, which is close to the initial value —5.12% of the “shadow rate” 
found in (Gorovoi and Linetsky, 2004) for Model-1. However, due to the discrete 
spectrum in Model-2, the fact that r_ is negative is not particulary informative. 
The short-rate modes that contribute into the bond (and other) prices are given 
by Xn = ^0 + En (recall that in this case xif) = ^o, so Xn(t) are constant), which 
are positive in this model. The hrst ten values of Xn are given in column 1, Table 
[3l It is clear that at long maturities only a few lowest-lying levels have signihcant 
contributions int o fl56ll . However, at short maturities a significant number of levels 
must be includedc/j 

The Model-2 and empirical yields are plotted in Figure 1. The fit in Model-2, 
which has fewer (to wit, 3) parameters than Model- lo is actually better than in 

Equivalently, we can simply set < = 0, so x(t) = tq. 

Typically, there are not that many maturities available, so reconstructing x(s) = dsr]{t, s) + 
x(t) and i^(s) = x(s) would involve piecewise polynomial splines. 

The Airy functions Ai{y) and Bi{y) and their zeroes are built-in within the “gsl” package in 
R. Their integrals are not, so we evaluated them (see (1751) . Appendix lAT) . 

A similar model was independently discussed by Rogers (1995). 

All dimensionful parameters are quoted in the units of 1 year. Note that r*, ro, z, /? and En 
have dimension (time)“^, a has dimension (time)”'^/^, and Wt has dimension (time)^/^. 

In our computation it was (more than) sufficient to truncate the series in (15611 at n = 300. 

In Model-1 the underlying “shadow rate” process Xt follows the Vasicek model ([5|) with 
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Model-1. The RMSE between the empirical yield curve (Tabled! column 4) and the 
Model-1 yield curve (Tabled! column 5) RMSE(Model-l) ~ 6.37 x 10“^. The RMSE 
between the empirical yield curve and the Model-2 yield curve (Table d! column 6) 
RMSE(Model-2) 5.91 x lO”"^. This implies that the mean-reverting feature in 
Model-1 (which introduces two additional parameters) apparently does not improve 
the £t. We have also computed a straightforward (not piecewise spline) cubic £t with 
the intercept, which amounts to fitting a general cubic polynomial into the empirical 
yield curve using a linear model. The RMSE between the empirical yield curve and 
the cubic fit yield curve (Table d! column 7) RMSE(Cubic) ~ 6.60 x 10“^. From 
Tabled! column 7 it is evident that both Model-1 and Model-2 are significantly better 
than the cubic fit. The cubic fit works very well at long maturities, but fails badly 
at short maturities. This is not surprising considering that, as mentioned above, at 
long maturities only the lowest few modes contribute, whereas at short maturities 
a large number of modes do. Furthermore, this suggests that a good model fit 
(both for Model-1 and Model-2) may not be a universal feature and may not persist 
to other data, in which case the “residual” yield might have to be attributed to 
nontrivial drift. The “residual” drift in Model-2 is plotted in Figure 2. 

3.3.3 An Illustrative Example: Recent US Treasury Yield Curve 

Table |21 summarizes the Model-2 fit for a recent US Treasury yield curve. The 
calibration procedure is the same as above. Column 3 of Table |2! corresponds to 
the Model-2 fit based on all maturities T; the empirical and Model-2 yields are 
plotted in Figure 3, and the “residual” yield is plotted in Figure 4; the first 10 
values of the Xn spectrum are given in Table [3! column 2. Column 4 of Table [21 
corresponds to the Model-2 fit based only on the maturities T > 1 yr; the empirical 
and Model-2 yields are plotted in Figure 5, and the “residual” yield is plotted in 
Figure 6; the hrst 10 values of the Xn spectrum are given in Table [3! column 3. For 
comparison purposes, columns 5 and 6 of Table [2! contain straightforward cubic hts 
into the empirical data for all maturities and T > 1 yr maturities only, respectively. 
We have the following RMSE for the above hts: RMSE(Model-2, All T) 1.99 x 
10-3; RMSE(Model-2, T > 1 yr) ^ 4.91 x 10"^ RMSE(Cubic, All T) ^ 5.07 x 10“^ 
RMSE(Cubic, T > 1 yr) 5.39 x 10“^. The Model-2 fit, which assumes vanishing 
drift, is better than the cubic fit for longer maturities, but not for shorter ones. 

Looking at the results it is clear that Model-2 provides a very good fit for maturi¬ 
ties T > 1 yr, but not for short maturities T < 1 yr. For these short maturities, the 
“residual” yield attributed to nontrivial drift is large and cannot be neglected. More¬ 
over, apparently, the drift is not even approximately constant - for a constant drift 
i>{t) = uq we would have the “residual” yield of the form Rr{t,T) = (z^o/2) {T — t). 
So, the drift appears to have nontrivial time dependence. 

constant tr, 9 and a, which together with z = rt give 4 parameters. The fifth parameter is the 
location of the sewing point r_, i.e., n = + (W)'*". In (Gorovoi and Linetsky, 2004) it is set to 

zero, r_ = 0. However, here too the spectrum can be nonnegative for a range of r_ < 0. 


12 



3.4 Bond Pricing on an Interval 

For the sake of completeness, let us briefly discuss zero-coupon T-Bond pricing in 
the Ho and Lee model on an interval. We have 


TT 


Cn= dx'ipnix) = - an[(j)n{0) Gi'{-en) - 4>n{L) Gi'{a L - Cn)] (65) 
.In O' 


where we have used fl7T]) and fl72D (see Appendix |^, and (j)n{x) is defined in fH7)l . 
The T-bond pricing function therefore is given by the following simple formula: 


v{z,t,T) = (66) 

n=l \ ^ J 

, _ 0n(O) Gi'{-en) - 0n(L) Gi' {a L - e„) 

e„ 02 (o) + («L-e„) 02 (L) ^ ^ 

where Xnif) = xif) + En = xif) + ^n, and e„ > 0 are the roots of fH9|) . Recall that 

a = (2ct)^/^ and /? = (ct^/ 2)^/^. If we assume vanishing drift, then the Ho and Lee 
model on an interval has four parameters to fit: z, a, tq (the lower reflecting barrier 
is at r_ = tq) and L (the upper reflecting barrier is at r_|_ = rg -|-T). The calibration 
can be done as above with the caveat that numerical estimations involving the Airy 
function Bi{y) are trickier as Bi{y) diverges for large positive y (see Appendix lAI 
for some useful formulas). 


4 Concluding Remarks 

The condition fl 2 ^ restricts the allowed drifts iz{t). As mentioned above, the l.h.s. 
of this condition should be some positive number, which we will denote by Rmin- 
Then we havj^ 

X*if) + xif) + El > Rmin ( 68 ) 

Let us now consider constant drift ^{t) = v. We have x(f) = ro -|- z/ t, so ri{t,T) = 
z/ (T — t)^ /2 and x*{t) is +oo if z/ > 0 and —oo if z/ < 0. This implies that constant 
negative drift is not allowedly This is evident from the fact that when the drift 
is constant we have = ViWt) + v t. In fact, positive constant v is not realistic 
either because it would drive rt to large values. A reasonable assumption then is 
that x*if) should be finite, which implies that 

hm|x(s)|<x (69) 

S—)-CXD 


where y is finite. 

This holds for general rt = V{Wt) + x(i). If V{x) is confining {e.g., V{x) = 2:^/2), then 

there is no need for a reflecting boundary as the spectrum is discrete and bounded from below. 

In the solution corresponding to introducing a reflecting barrier for rt (as opposed to 
Wt) constant negative drift a priori is allowed as r* is reflected when it hits its lower boundary r*. 
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Another point worth commenting on is that, while for time-homogeneous cases 
{i.e., for time-independent parameters) we can consider “option-like” models with 
Tt = and the underlying “shadow” rate Xt following processes such as Vasicek, 

CIR, etc., which have discrete spectra, Xt cannot follow the Ho and Lee model as 
the latter has continuous spectrumQ In contrast, we can have a reflecting barrier in 
the Ho and Lee model. Such a barrier can be introduced for Wt for a time-dependent 
drift (subject to the above conditions), and also on for constant drift, which is 
negative in the case with only the lower barrier, but can be positive in the case 
where both lower and upper barriers are present. However, constant drift is limited 
in its applicability as can be seen from the U.S. Treasury data we discussed above. 

A Some Properties of Airy and Scorer Functions 


In this appendix we collect some properties of the Airy and Scorer functions used in 
the main text. Unless stated otherwise, these properties are taken from (Abramowitz 
and Stegun, 1964). 


• Integral identity: 


Jdy fi{y) f 2 {y) = y fi{y) f 2 {y) - fi{y) 

(70) 

where each /i and /2 can be either Ai{y) or Bi{y). Eq. flTUD can 
integration by parts and using the Airy equation f”{y) = y fi{y), i 
• Integral identities: 

be obtained via 
= 1,2. 

1 

— dw Ai{w) = Ai{y) Gi'{y) — Ai'{y) Gi{y) 

^ dy 

(71) 

— [ dw Bi{w) = Bi'{y) Gi{y) — Bi{y) Gi'{y) 
TT Jo 

(72) 

fOO ^ 

/ dw Ai{w) = - 
Jo 3 

(73) 

where Gi{y) is a Scorer function. Let 


A{—y) = f dw Ai{w) 

J-y 

(74) 

B{—y) = 1 dw Bi{w) 

J-y 

(75) 

where ?/ > 0. The integrals A{—y) and B{—y) can be evaluated numerically. It 
follows from fl73D that 

j dw A{w) = - + A{-y) 

(76) 


Here one can wonder if it would make sense to consider “option-like” models for W*. However, 
if we replace Wt by (Md)+ in rt = V{Wt) + x(t), the spectrum will be continuous. 
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( 77 ) 


where |/ > 0. It further follows from fl7T]) and d72|) that 

^ Ai(-p) At'(-p) 

which allows to evaluate Gi'{—y) numerically for |/ > 0. In some cases, the integral 
representation for Gi{y > 0) given in (Gil, Segura and Temme, 2001) may be useful 
for evaluating Gi'{y > 0). 

• Leading asymptotic behavior 

Ai{-y) ~ sin Q y^^^ + (78) 

Gi'{-y) ~ ^ sin Q y^^^ + (79) 

at 1. 

• Leading asymptotic behavior at large n: 

/ Ovr \ 

('—1 , 

(si) 

\Ar 

( 82 ) 

Vvr 

where (0 > > ^2 > • • •) are the zeros of Ai\y). 

B Reflecting Barriers for the Short-rate Process 

In this section we discuss the issues associated with introducing reflecting bound¬ 
aries directly for the short-rate process as opposed to the underlying Brownian 
motion process Wt- Our starting point is ([5]) together with the terminal condition 
v{z,T,T) = Y(z). When a{rt,t) and p{rtA) have no explicit time dependence, re¬ 
flecting barriers for r* can be introduced as in, e.g., (Goldstein and Keirstead, 1997) 
and (Gorovoi and Linetsky, 2004) @ Even if a{rt, t) has no explicit time dependence, 
for general, explicitly time-dependent ^{rt, t) things are trickier. If fact, for our pur¬ 
poses here, to illustrate our main point, it will suffice to consider the case where 
a{rt, t) = a is constant and z/(rt, t) = vit) is independent of rt, but is a deterministic 
function of t. Eq. (jS]) then reduces to 

z/(t) dzv{z, t, T) + dtv(z, t, T) -|- -cr^ d‘^v{z, t,T) — z v{z, t,T) = 0 (83) 

In this case, (I5|) can be transformed into the Schrodinger equation with time-independent 
potential (a.k.a. Liouville normal form) via a Liouville transformation. 
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An analogue of the Liouville transformation in this case is 


v{z,t,T) = exp[-^^^^] u{z,t,T) 




Then the PDE for u{z,t,T) reads: 
dtu{z,t,T) + dlu{z,t,T) - 


z 1 + 


2a2 




+ 


( 84 ) 


u{z, t,T) = 0 (85) 


This is the Schrodinger equation with time-dependent potential - unless z/(t) is 
constant, that is. When z/(t) is a linear function of t so that i>{t) = p is constant, 
the potential is time-dependent but factorized, so in this special case the solution 
to is given by 

u{z,t,T) = exp ds w{z,t,T) (86) 

where w{z,t,T) solves the following Schrodinger equation with a linear potential: 

dtw{z, t, T) dlw{z, t,T) - z (^1 + w{z, f, T) = 0 (87) 

which can be solved as in the main text; provided, however, that the boundary 
condition can be consistently specihed. 

This is where a difficulty arises irrespective of whether v{t) is a linear function 
of t or not - so long as it is not constant, that is. Thus, if we impose a reflecting 
(Neumann) boundary condition for some = r* 


d,v{z,t,T)\^^^^ = 0 

then for u{z,t,T) we have the following Robin boundary condition: 


dzu{z,t,T)\^^^^ = 


u{t) 




u(z,t,T) 


( 88 ) 


( 89 ) 


For constant v{t) we can perform the separation of variables. However, even for 
linear v{t) the straightforward separation of variables procedure cannot be applied 
and the problem becomes more difficult to solve. 

Let us briefly outline the solution for constant u{t) = u, primarily to illustrate 
the difference with the corresponding solution in the case of a reflecting boundary 
for the underlying Brownian motion we discuss in the main text. We havj^ 


v{z,t,T) ='il:{z//3,T-t) (90) 

CX) 

'ip{y,T) = Ai{y - Cn) (91) 

n=\ 

Here we focus on the solution with only one (lower) boundary. The solution with both (lower 
and upper) boundaries can also be readily constructed and involves both Ai{y) and Bi{y). 
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where A„ 
equation 


(3 ( 7 ^ + e„), (3 = (ct^/2 )^/^, 7 = z//(2ct^)^/^, and e„ are the roots of the 
Ai' (?/* - e„) = 7 (?/* - e„) (92) 


which is a consequence of the boundary condition fl89|) corresponding to a reflecting 
boundarj@ at = r*, and = r*//5. The coefficients dn are hxed via the terminal 
condition v{z,T,T) = Y(z). For a zero-coupon T-bond we have Y{z) = 1 and 


^ ^ j^^dye^yAi{y-en) ^ /(7, -?/.) +/(q) 

47 dy [Aliy - e„)]^ (e„ - y^ - 7 ^) - e^)]^ 

I('y,y)= f dw e^^Ai(w) 

J-y 

^ ^00 

1(7) = / dm e^'^Ai(w) = 

Jo 





34/3 r(4/3) 




2 . 5._T\^2 
3’ 3’ 3 M 


3^5/3 r(5/3) 


(93) 

(94) 


(95) 


where r(|/) is the Gamma function, and iTi(a; b]y) is a generalized hypergeometric 
function. In (195]) we have used a Laplace transform of the Airy function Ai{y), see 
Eq. (9.10.14) in (DLMF, 2015). In practice, if I 7 I is not large, then the integral 
1 ( 7 ) can be evaluated by truncating the series 


7(y) = iV 

3^F(n/3 + l) 


(96) 


which follows from the Mellin transform of Ai{y), see Eq. (9.10.17) in (DLMF, 2015). 
Asymptotically, this series behaves as the Taylor expansion of (1/3) exp( 7 ^/ 3 ). 
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Table 1: The calibration results for the zero-coupon T-bond pricing formula flSB]l 
in the Ho and Lee model on a semi-inhnite line with vanishing drift (“Model-2”). 
The hrst hve columns are taken from (Gorovoi and Linetsky, 2004). The hrst three 
columns give the Bloomberg data as of t = 02/03/2002 for the Japanese Government 
Bonds: coupon, maturity T and price. The fourth column gives the corresponding 
bootstrapped zero-coupon bond yields. The hfth column gives the calibrated model 
yields calculated in (Gorovoi and Linetsky, 2004) for Black’s model of interest rates 
as options (Black, 1995) with the underlying “shadow rate” following the Vasicek 
model (“Model-1”). The sixth column gives calibrated model yields we calculated 
based fl56|l . The calibrated model parameters are 2 ; ~ —0.00184, (3 0.0924 (which 

implies a ~ 0.0397), tq ^ —0.05834. For illustrative purposes we have kept 3 digits 
in the sixth column. In the seventh column we have included the yield curve from 
htting a general cubic polynomial into the empirical data in the fourth column. 


Goupon 

Maturity 

Price 

Empirical 
0-Goupon 
Yield (%) 

Mo del-1 
Yield (%) 

Model-2 
Yield (%) 

Gubic Fit 
Yield (%) 

4.2 

3/20/2003 

104.648 

0.02 

0.03 

0.023 

-0.106 

3.4 

3/22/2004 

106.900 

0.14 

0.17 

0.106 

0.156 

4.4 

3/21/2005 

112.729 

0.30 

0.36 

0.338 

0.395 

3.1 

3/20/2006 

110.481 

0.54 

0.57 

0.571 

0.615 

2.6 

3/20/2007 

109.326 

0.76 

0.78 

0.788 

0.818 

1.9 

3/20/2008 

105.578 

0.98 

0.98 

0.988 

1.005 

1.9 

3/20/2009 

104.723 

1.24 

1.16 

1.169 

1.175 

1.7 

3/22/2010 

102.521 

1.40 

1.33 

1.333 

1.331 

1.4 

3/21/2011 

99.314 

1.51 

1.48 

1.481 

1.472 

1.5 

12/20/2011 

99.997 

1.53 

1.59 

1.584 

1.569 

3.8 

9/20/2016 

123.287 

2.11 

2.09 

2.084 

2.040 

2.1 

12/20/2021 

98.411 

2.29 

2.44 

2.434 

2.360 

2.4 

11/20/2031 

94.810 

2.88 

2.79 

2.801 

2.869 
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Table 2: The calibration results for the zero-coupon T-bond pricing formula fl56|l in 
the Ho and Lee model on a semi-inhnite line with vanishing drift (“Model-2”). The 
hrst two columns are taken from the U.S. Treasury website (U.S. Treasury, 2015) as 
oft = 01/29/2015: maturity T and empirical yield Re{t, T). The third column is the 
Model-2 fit including all maturities, for which the calibrated model parameters are 
2 : ~ —0.0027, f5 ~ 0.2516 (which implies a ~ 0.1785), tq ^ —0.23163. The fourth 
column is the Model-2 £t including only the maturities of 1 year and longer, for 
which the calibrated model parameters are z ~ 0.0012, (3 0.2085 (which implies 

a ~ 0.1346), ro ~ —0.1879. For illustrative purposes we have kept 3 digits in the 
third and fourth columns. In the hfth (all maturities) and sixths (only the maturities 
of 1 year and longer) columns we have included the yield curve from htting a general 
cubic polynomial into the empirical data in the second column. 


Maturity 

T 

Empirical 
Yield (%) 

Model-2 
Yield (%) 
(All T) 

Model-2 
Yield (%) 
(T > 1 yr) 

Cubic Fit 
Yield (%) 
(All T) 

Cubic Fit 
Yield (%) 

(T > 1 yr) 

1 mo 

0.01 

0.046 

— 

-0.040 

— 

3 mo 

0.03 

0.524 

— 

0.016 

— 

6 mo 

0.07 

-0.005 

— 

0.097 

— 

1 yr 

0.17 

-0.151 

0.209 

0.252 

0.224 

2 yrs 

0.51 

0.292 

0.464 

0.539 

0.521 

3 yrs 

0.84 

0.733 

0.797 

0.795 

0.785 

5 yrs 

1.28 

1.314 

1.295 

1.221 

1.224 

7 yrs 

1.59 

1.629 

1.594 

1.544 

1.554 

10 yrs 

1.77 

1.880 

1.845 

1.865 

1.878 

20 yrs 

2.11 

2.175 

2.147 

2.086 

2.078 

30 yrs 

2.33 

2.273 

2.249 

2.335 

2.336 
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Table 3: The first 10 values Xn contributing to the zero-coupon T-bond pricing 
formula fl56|) in the Ho and Lee model on a semi-inhnite line with vanishing drift 
(“Model-2”). The hrst column corresponds to the Japanese Government Bond data 
in Table [U (JGB). The second column corresponds to the US Treasury yield data 
for the £t including all maturities in column 3 of Table[2] (USTl). The third column 
corresponds to the US Treasury yield data for the fit including only the maturities 
of 1 year and longer in column 4 of Table |2] (UST2). Note that Xn are independent 
of t when the drift vanishes. 


Xn 

JGB (%) 

USTl (%) 

UST2 (%) 

Xi 

3.578 

2.470 

2.451 

X2 

24.187 

58.562 

48.934 

X3 

38.718 

98.111 

81.708 

X4 

51.134 

131.906 

109.714 

X5 

62.309 

162.321 

134.919 

X6 

72.628 

190.407 

158.194 

X7 

82.306 

216.749 

180.023 

X8 

91.478 

241.713 

200.711 

X9 

100.236 

265.549 

220.464 

Xio 

108.646 

288.438 

239.432 
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Figure 1. Empirical (circles) and Model-2 (triangles) yield curves corresponding to the 
data in Table [H 
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Time to maturity (years) 


Figure 2. The “residual” yield attributed to nontrivial drift in Model-2 corresponding to 
the data in Table [TJ 
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Figure 3. Empirical (circles) and Model-2 (triangles) yield curves corresponding to the 
data in Table [2] for all maturities. 
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Time to maturity (years) 


Figure 4. The “residual” yield attributed to nontrivial drift in Model-2 corresponding to 
the data in Table [2] for all maturities. 
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Figure 5. Empirical (circles) and Model-2 (triangles) yield curves corresponding to the 
data in Table [2] for the maturities of 1 year and longer. 




Time to maturity (years) 


Figure 6. The “residual” yield attributed to nontrivial drift in Model-2 corresponding to 
data in Table [2] for the maturities of 1 year and longer. 
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